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We present a new technique to extract the cosmological information from high-redshift supernova 
data in the presence of calibration errors and extinction due to dust. While in the traditional 
technique the distance modulus of each supernova is determined separately, in our approach we de- 
termine all distance moduli at once, in a process that achieves a significant degree of self-calibration. 
The result is a much reduced sensitivity of the cosmological parameters to the calibration uncer- 
tainties. As an example, for a strawman mission similar to that outlined in the SNAP satellite 
proposal, the increased precision obtained with the new approach is roughly equivalent to a factor 
of five decrease in the calibration uncertainty. 
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I. INTRODUCTION 

The study of type-la supernovae provided the first indication that the expansion of the Universe is accelerating 
This conclusion is supported by the combination of Cosmic Microwave Background results ^] and large-scale structure 
measurements Results from more recent supernova surveys further strengthen the evidence for the accelerated 
expansion of the Universe and the necessity for some mysterious mechanism that drives it, code-named "dark energy." 

Supernova magnitudes at peak brightness and redshifts are the ingredients that allow their use as cosmological 
probes. Peak magnitudes can be corrected with the "stretch" 0], Amis 0, or MLCS techniques to turn type-la 
supernovae into standardized candles to a level around 0.10-0.15 mag (10-15% in flux). Therefore, measuring the 
supernova apparent magnitude provides the distance and hence the lookback time to when the explosion occurred, 
while measuring its redshift z provides the scale factor a = {1 + z)~^ at the time the supernova light was emitted. 
These two statistics for each supernova from a set that spans a range of redshifts provide a mapping of the cosmic 
expansion with time. 

Neglecting for the moment any secondary effects, one can relate the apparent peak magnitude m(z) of a supernova 
at redshift z to the luminosity distance d{z) = r{z){l + z), where r{z) is the dimensionless comoving distance, through: 



m(z) = -2.5 logio (T/To) = M+ /i(z) =M + 5 log^o d{z) + 25 



(1) 



where J- is the observed supernova flux and J-q is the flux corresponding to zero magnitude of the magnitude system. 
Ai = M — SlogiQ [TJo/lOOkm s~^Mpc~^] , can be regarded as an unknown nuisance constant, M being the absolute 
supernova magnitude, while /z is the dimensionless supernova distance modulus. With this definition r(z) is also 
dimensionless and does not depend on Ho: 
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where we have assumed a fiat universe and w{z) is the equation of state of the dark energy, w{z) = —1 for a 
cosmological constant. In more general models w{z) would simply be a function that captures the effect of dark 
energy. Following llOj we consider the parameterization 



w{z) = Wo + Wa {I - a{z)) 



(3) 



where wq and Wa are constants that determine, respectively, the dark-energy equation of state now and its evolution. 
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In a real observation there are many more effects that have to be accounted for: weak lensing, Milky Way extinction, 
etc. In this paper, we will concentrate on two (inter-related) effects: host-galaxy extinction due to dust and mis- 
calibration. Some of the other effects are treated in . Dust in the supernova host galaxy will dim the light coming 
from the supernova, therefore biasing its measured magnitude toward larger values. However, ordinary dust dims 
light of different wavelengths differently, absorbing more blue light than red light. By measuring the supernova flux 
in several broadband filters, one can determine the reddening and hence the overall dimming that dust creates. In 
order to proceed with the last step two things are needed: a good understanding of the gross features of non-extincted 
type-la supernova spectra, and a knowledge of the relationship between differential dimming across wavelengths and 
absolute extinction. We will assume a p erfect knowledge of the former and for the latter we will use the standard 
Cardelli, Clayton and Mathis model ^!z|. In this model, the extinction A (that is, the dimming in magnitude units) 
at a wavelength A can be written as: 

A{X) = Av ■ (a(A) + h{\)/Rv) (4) 

where the parameters Ay and Rv control, respectively, the overall amount of dimming and its dependence on wave- 
length, and a(A) and 6(A) are known functions. The value of Ay will be different for different supernovae depending 
on the properties of the host galaxy and the position of the supernova inside it. The value of Ry will depend on 
the type of dust (grain size, for instance), and observations show |T3 that it varies significantly across redshift and 
between galaxies. In the following, we will assume no prior knowledge on either Ay or Ry. Then, given enough 
independent measurements of the magnitudes of a given supernova i in several broadband filters centered around 
wavelengths Aq, one can in principle determine directly from the observations the dust parameters Ay, and Ry for 
supernova i, together with its distance modulus. 

Filter calibration uncertainties affect the measured magnitude in each passband. Therefore, they affect not only the 
overall peak magnitude but also the extinction correction that, as mentioned above, depends on understanding the 
supernova colors (defined as the ratio of the supernova flux in two different filters). The calibration error in a given 
observer filter will affect each supernova in a different (but correlated) way, since a given passband in the supernova 
rest frame will be observed with a different filter in the observatory depending on the supernova redshift. Calibration 
uncertainties are going to be among the limiting uncertainties lllj in the new generation of supernova surveys being 
currently planned. 

In the standard way of extracting cosmological information from supernova data, the data for a given supernova i 
is used to determine the values of /i*, Ay, and Ry and their uncertainties. Then the calibration error is propagated 
through that supernova's measurements. As a result, a non-diagonal covariance matrix for the set of distance moduli 
{/i*} is obtained. The measurements of all the {/i*} and their covariance matrix are then used in a fit to determine 
the cosmological parameters fi™, wq, and Wa, using Eqs. (0) and 

In our approach, all {/i*} are determined simultaneously. In this way, all information available is used in an 
optimal way, which achieves a certain degree of self-calibration. The cosmological parameter fit then follows as 
above. Equivalently, one can just bypass completely the step of extracting the distance moduli and determine directly 
the cosmological parameters from the set of magnitudes measured for all supernovae {m* } (« runs over supernovae, 
whereas j runs over passbands in the supernova rest frame). In Section^ we present a fiducial supernova mission: 
the dataset it provides and the algorithms through which the dataset is used to measure the dark energy. Within that 
framework, we elaborate further on the calibration error and its interplay with extinction. Section [llll describes both 
the traditional data analysis method and our new approach. Results using both approaches are given in Section ll VI 
for a fiducial mission and the statistical advantage that the new, optimal method affords becomes apparent. Finally, 
we draw our conclusions in Section Ivl 

II. CALIBRATION AND SUPERNOVA COSMOLOGY 

Although our results have a more general validity, in the remainder we discuss a specific strawman mission along 
the lines of the SNAP satellite proposal 0| . The mission observes 2000 type- la supernovae over a redshift range from 
up to 1.7 with a set of nine logarithmically-spaced broadband filters centered at wavelengths Aq, = Aq • (1 -I- 0.16)", 
a = 0, . . . , 8, with Aq = 440 nm corresponding to the B band. For every supernova, the analysis is based on the B 
passband in the supernova rest frame and any other band redder than this which, after redshifting to the observer 
frame, still ends up overlapping one of the nine filters. So only optical and infrared information in the supernova rest 
frame is used. The current lack of understanding of supernova heterogeneity in the UV prevents the robust usage 
of shorter wavelengths (though see [isl Il6| for perspectives for using UV light curves to aid in measuring supernova 
distances) . 

Therefore, 9 — fc filters are used for a supernova with z between and z^+i, where Zk = (1 + 0.16)^^ — 1, k = 0, . . . , 6. 
With this choice of filters, at least three bands are available for each supernova. The calibration uncertainty enters 
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through the zero points of the filters, = 2.5 log^Q ^ao, a = 0, ... ,8, where TaQ is the observed flux measured in 
filter a for a standard object of known magnitude 0. So in order to specify a calibration error, one has to specify the 
9x9 covariance matrix for {-Z^}. 

We are now in a position to write the equation that will relate the measured light-curve-shape-corrected P, Q, 
and K-corrected | l7j| peak magnitudes m* with the parameters of our model (details on K-correction uncertainties, 
which we will not consider here, can be found in flS]): 

™}+2a(»j) =A^j+/i''(2';f^m,wo,Wa)-f AV-a(Aj)-|-Bi,-6(Aj) . (5) 

Indices i and j run respectively over supernovae and passbands. M.^ is M for passband j: M.j = Mj — 
5 logj^Q [i/o/lOOkm s^^Mpc , where Mj is the absolute magnitude of a type-la supernova measured in passband j. 
By is defined as Ay jRy and makes the problem linear on the dust parameters. Finally, Q;(i,j) is the index of the 
filter in which, after K-correction, rest frame passband j of supernova i is measured. For Zfc < < -Zfc+i, we have 
a{i, j) = j + k, and the index j runs from to 9 — — 1, so that there are 9 — k measurements for supernova i. 

The left-hand side of Eq. Q contains the measurements, while in the right-hand side we have the parameters to be 
determined: Aij, j = 0, . . . ,8, and /i'. Ay, By for each supernova. Alternatively, the fit parameters /i* can be replaced 
with the fit parameters f2m, wq, and Wa with /z' — s- ^m, wq, Wa) where the latter is the predicted magnitude as a 
function of the cosmological parameters. A careful examination of Eq. ((S} reveals that there are two non-independent 
parameters. Specifically, a simultaneous shift in all Ay to Ay + S, can be compensated by a shift of all A4j to 
Mj — S ■ a{Xj). Analogously, a change By By -I- 77 is exactly compensated by the shift Mj Mj — rj ■ b{Xj). 
Therefore, one has to eliminate two parameters. To do so, we assume that for a nearby supernova we know that it has 
no dust (for instance because it lies in an elliptical galaxy), or equivalently that its dust extinction is perfectly known. 
For our analysis, we choose to assume that for the supernova with index i = it holds Ay = By = 0. Furthermore, 
for simplicity, we assume that we know the peak magnitudes for that supernova, m^, perfectly well. Results do not 
depend strongly on this assumption, which can be justified because the statistical errors in the photometry of a nearby 
supernova are going to be much smaller than those for the high-redshift ones. We can now subtract from Eq. (jSJ the 
particular equation for i = 0, to get 

m} - m° + - Zj = fi'{z'; wq, i«a) - / + A\, ■ a(A^) + B'y ■ b{\j) , (6) 

where we have taken into account Za^Q j) — Zj when supernova has a low redshift < zi — 0.16. While the 
m'j magnitudes are assumed known, fjP is not and it becomes an additional free parameter. Now the index i runs 
from one to the number of supernovae minus one. The number of free parameters has been decreased by two with 
respect to Eq. Written in this form, it also becomes clear that only eight color calibrations {Za — Zp) rather than 
nine absolute calibrations are necessary for the cosmology determination. This is the framework that we use in the 
following. 



III. SIMULATED DATA ANALYSIS 
A. Supernova-by-supernova procedure 



In standard supernova-cosmology analysis, the distance to each supernova is determined independently. In our 
implementation of the traditional analysis procedure, values for Ay, By, and /i* — are determined from the set of 
magnitudes {'ti*} for supernova i. Only statistical errors crj g^^^^ in the measurements of the magnitudes are included in 
the fit. A Monte Carlo process is then used to compute the additional covariance matrix for the set of {/x'} due to the 
zero-point calibration uncertainties a a. cai, where the index a labels the filters. Finally, the cosmology fit is performed 
on all {/i*} using the covariance matrix that includes calibration errors plus an additional intrinsic dispersion (7?^^ 
for each supernova distance modulus that takes into account known supernova-to-supernova variability. The free 
parameters in the cosmology fit are , Vl^m w^, and Wa- For the fit we use the standard package minuit ^ . We have 
checked that a Fisher matrix calculation provides compatible error estimates. 



|http : //wwwlnf o . cern ■ ch/asdoc/minult/minmain ■ htm~ 
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B. Simultaneous procedure 

In the new approach, aU supernova data are used at once in order to determine all Ay, By, and /i* — /x". The 
input covariance matrix for the to* includes statistical errors (which are diagonal), calibration error and intrinsic 
dispersion. For a large simulated data set of about 2000 supernovae, this treatment is rather impractical, since the 
input covariance matrix is about 12000 x 12000 and non-diagonal, which makes it difficult to store and invert efficiently. 
Instead we can treat both the intrinsic dispersion and the calibration error as new free fit parameters, constrained 
within the known uncertainties. Then one can build a function whose covariance matrix is diagonal: 



TO,,- 



^ j ""j ~ "^j ~ "^3 ~ ^ ' ^V'"^ ' ^aj j + ^ I Ji_ j + ^ ( ^" j (7) 

fi' {z'; Qm, WO, Wa) - M° + • a i^j) + B'y ■ b + S'^ + Zj , (8) 



where in Eq. Q, again the indices i and j run over supernovae and passbands respectively, a runs over filters, and 
S'* is the new free parameter that accounts for the intrinsic supernova dispersion, one new parameter per supernova. 
Now the {Za} are considered as nine additional free parameters, constrained to vary within the calibration precision. 
The uncertainties cr* al^^, and aa, ced are as defined in Section llll Al Typically, for 2000 supernovae, we now have 
about 12000 measurements and over 8000 unknowns. After all {/i* — /i'^} are determined, the cosmology fit proceeds 
as in Section [ill Al but now using the {/i'} covariance matrix from the overall fit to all distance moduli. 

Alternatively, and equivalently, one can do without the determination of the distance moduli, and just fit all the 
magnitudes directly to the cosmological parameters and the other nuisance parameters, by using again the in 
Eqs. (|7I8|) . Now the free parameters in the fit are all the Ay, By, and Z^, as well as Sim, wq, and Wa- This 
is the final approach that we have used. A Fisher matrix calculation following these lines runs on a 64-bit^, 3 GHz 
Intel Xeon processor in about one hour for 2000 supernovae. 



IV. RESULTS 



A. Calibration models 



In this section we are going to present results based on two extreme calibration models. We will not attempt to 
construct a realistic model from observations of standard stars or calibrators. Instead, we will specify the calibration 
model by its 9x9 zero-point covariance matrix. 

In the first model, we assume that the zero-point covariance matrix is fully diagonal. While this is surely an 
unrealistic scenario (the same fundamental calibrator is probably going to be used for more than one filter, if not 
all), there will still be a diagonal component in the full calibration covariance matrix, originating for instance from 
Poisson uncertainties in the measurement of the fundamental calibrator. For simplicity, we will take the calibration 
errors in all filters to be equal, so that the zero-point covariance matrix will be simply Va/3 = cr^^i ' ^ap- We will call 
this model the "Diagonal" model. 

The other model we consider goes to the opposite extreme, and assumes that all sources of error during the 
calibration process come from the uncertainty in a single parameter, such as the temperature of a hot white dwarf 
used as a calibrator. In this case the zero-point parameters are replaced with a single temperature parameter T and 
Eq. O has to be modified to replace the last term with a term (T — Tq)'^ /a^ where Tq is the given temperature 
and ctt the precision with which it is known. Then, in Eq. ^ the zero points Z^ are turned into functions of the 
temperature parameter by integrating in the corresponding filter the black body spectrum of the white dwarf with 
that temperature. We will call this model the "Temperature" model. 

As mentioned above, these are two very extreme toy models for the calibration covariance matrix. The true 
covariance matrix will be much more complicated and will depend on the specific details of the mission and on the set 
of standard calibrators used. However, these two models should span the range of possible calibration error models. 



^ Use of 128-bit (double) precision arithmetic is crucial in the process of inverting the about 6000x6000 heavily-correlated Fisher matrix. 
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TABLE I: The redshift distribution N(z) of the supernovae employed in the fiducial survey. The redshifts z given in the table 
correspond to the value for all supernovae in that bin. 



z 


0.05 


0.17 


0.35 


0.57 


0.82 


1.11 


1.50 


Niz) 


317 


82 


219 


412 


441 


427 


400 



TABLE II: Uncertainties in the cosmological parameter determination in the two analysis methods for several values of the 
calibration uncertainty in the "Diagonal" error model. 



0"cal 




SN by SN 


Simultaneous 


0.000 


a{wo) 


0.064 


0.064 




cr{wa) 


0.30 


0.30 


0.001 


a{wo) 


0.082 


0.068 




cr{wa) 


0.40 


0.33 


0.005 


a{wo) 


0.099 


0.071 




cr{wa) 


0.59 


0.43 


0.010 


a (wo) 


0.11 


0.075 






0.81 


0.53 



B. No calibration error 

To start with, we consider the case with no calibration error. That is, we fix all parameters {Z^} in Eq. Q or 
equivalently, we take the limit (Ja, ced 0. Let us now specify completely our fiducial survey. Table^gives the number 
of supernovae in each redshift bin, based on the numbers in llj. All supernovae in a bin have been assigned the 
same actual redshift, which is also reported in the table. In total there are 2298 supernovae, 300 of which would be 
obtained at low redshift by ground-based observations like those of the Supernova Factory [13 • We consider that the 
calibration error in the observation of these 300 additional supernovae is completely correlated with the calibration 
error in the original 2000 supernovae. This is what will happen if the dominant calibration uncertainty comes from 
the limited knowledge of the intrinsic properties of the calibrators, and the same calibrators are used in both surveys. 

For the statistical errors in the magnitudes we choose cr^^ gtat = 0.01 for all passbands in all supernovae. This 
uncertainty encompasses uncorrelated statistical uncertainties introduced from the photometry process (e.g. flatfield 
errors, Poisson noise) and intrinsic uncorrelated supernova color dispersion. The intrinsic dispersion in the supernova 
magnitude is ctjj^^ = 0.15 for all supernovae, independent from supernova to supernova and fully correlated for all 
passbands within one supernova. Other than the supernova measurements, we will include in our synthetic data 
sample a measurement of the distance to the surface of last scattering from the future Planck survey, with a relative 
precision of 0.7%. Finally, we take as our fiducial cosmology the flat ACDM model, wq — —1, Wq — 0, with ri„j = 0.28. 

The resulting uncertainties in the cosmological parameters wq and Wa are given in the first two rows of Table ^ 
When there are no calibration errors the results obtained with the two analysis methods fully agree, as they should, 
given that in this case there is no information shared between any two supernovae, and therefore it should not make 
a difference whether they are analyzed separately or all at once. 

C. Diagonal calibration model 

For the "Diagonal" calibration error model, we will choose Vaf3 = cr^ai ' ^af3, with cTcai varying between 0.001 and 
0.01. The results again are in Table HTl Since now, because of the calibration error, there is correlated information 
shared among the supernovae, there is a clear advantage in analyzing all supernovae at once, to the point that, for 
instance, the effect of a 0.005 calibration uncertainty using the simultaneous fit is reduced to approximately the effect 
of a 0.001 calibration uncertainty when considering each supernova individually. 

Inspection of the final uncertainties in Za after the Fisher matrix calculation reveals the amount of self-calibration 
achieved: for example when 0.010 is assumed as the prior error in each Z^^ the final (posterior) Z^ errors are reduced 
to values ranging between 0.004 and 0.008, with large positive correlations between neighboring filters. 
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D. One-parameter calibration model 

In the "Temperature" model, one single parameter, in this case the temperature of a single Tq = 20000K hot white 
dwarf, defines the whole calibration error matrix. In this case, the final precision in the cosmological parameters is 
rather insensitive to the calibration error as was found in For instance, assuming a temperature error such that 
it results in an uncertainty for the zero point of the first filter of ctq, cai = 0.01 (and totally correlated uncertainties in 
the other zeropoints) results in an inappreciable change in the error on wq and Wa- Even if we assume cto, cai — 0.10, 
the uncertainties increase only marginally from 0.064 to 0.066 (wq) and from 0.30 to 0.31 (wa)- These results hold 
when performing the analysis in either of the two methods. The reason is that with one single parameter defining 
the calibration error, self-calibration can occur already within a single supernova, provided there are at least four 
passbands available, which is the case for over 80% of our supernovae. 

V. CONCLUSIONS 

In summary, we have presented a new method of analysis of high-redshift supernova data in the presence of 
calibration uncertainties and dust extinction that performs the statistical analysis of all data at once in order to 
self-calibrate the zeropoints of the filter set. For a fiducial mission inspired on that specified in the SNAP proposal, a 
significant reduction of the sensitivity of the cosmological parameters to the calibration uncertainty is achieved with 
the new method, equivalent to a reduction of about a factor five in calibration uncertainty. The advantage of the 
simultaneous analysis is not only applicable to calibration uncertainties, but to all sources of uncertainty that are 
non-trivially correlated between supernovae. 

Calibration uncertainties are no longer treated as irreducible; the supernova data provide an avenue for self- 
calibration of magnitude zeropoints. Assuming that type-la supernovae are standardizable candles in our passbands, 
their colors as a function of redshift can be predicted. Measurements of supernova colors at different redshifts there- 
fore constrain the differences in zeropoints, {Z^ — Zp}. Some of the color information is spent in correcting for dust 
extinction of individual supernovae, however given enough supernovae observed in enough passbands, a significant 
improvement in calibration is still achieved. 

This paradigm for simultaneous supernova analysis gives rise to new questions that can be addressed in future 
studies. The interplay between the supernova sample, input calibration uncertainty, cosmology priors, photometry 
uncertainty, and intrinsic supernova magnitude and color dispersion needs to be explored. The effects of merging 
disparate supernova samples whose calibrations are not perfectly correlated can be studied. Equally important is an 
exploration of systematic uncertainties incurred if incorrect models for either dust-extinction or intrinsic supernova 
colors are used. The generation of a realistic covariance matrix that faithfully describes the calibration process would 
show where between the "Diagonal" and "Temperature" models realistic observations will lie. 

For large data sets like those in the new method presents some practical numerical problems. In this study, a 
Fisher matrix analysis has been performed, which provides reliable estimates of the attainable precision. In the future, 
a Monte Carlo approach should allow for a real fit to the data, as well as facilitate the inclusion of non-Gaussian 
priors in, for instance, the distributions of Ay and Ry. 
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